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We present a method for simulating the time evolution of one-dimensional correlated electron- 
phonon systems which combines the time-evolving block decimation algorithm with a dynamical 
optimization of the local basis. This approach can reduce the computational cost by orders of 
magnitude when boson fluctuations are large. The method is demonstrated on the nonequilibrium 
Holstein polaron by comparison with exact simulations in a limited functional space and on the 
scattering of an electronic wave packet by local phonon modes. Our study of the scattering problem 
reveals a rich physics including transient self-trapping and dissipation. 
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Phonon degrees of freedom play an important role 
in the nonequilibrium properties of correlated materi¬ 
als. In particular, time-resolved spectroscopy [1, 2], pho- 
toinduced phase transitions [3, 4], and transport through 
low-dimensional or molecular junctions [5-7] call for the¬ 
oretical investigations of the nonequilibrium dynamics of 
charge carriers coupled to lattice vibrations. Quantum 
lattice models such as the Holstein model [8] are often 
used to describe the low-energy properties of strongly cor¬ 
related electron-phonon (EP) systems. Analytical stud¬ 
ies of these systems out of equilibrium are very difficult 
and reliable results are scarce. Therefore, theorists of¬ 
ten turn to numerical methods to investigate them [9- 
11]. However, accurate numerical simulations of corre¬ 
lated lattice systems are very challenging because of the 
rapid increase of the Hilbert space dimension with system 
size and phonon number fluctuations. Similarly, comput¬ 
ing the nonequilibrium dynamics of correlated bosons is 
a significant challenge in a great variety of physical sys¬ 
tems such as nonlinear optical systems [12, 13], quantum 
dissipative systems [14, 15], and low-energy models of 
quantum chromodynamics [16]. 

In this paper we present a method for simulating the 
time evolution of one-dimensional (ID) lattice models 
with strongly fluctuating bosonic degrees of freedom for 
long periods of time. It combines the time-evolving block 
decimation (TEBD) [17, 18] with a local basis optimiza¬ 
tion (LBO) approach [19] to reduce the computational 
cost significantly. The key idea is to optimize the local 
bases for the bosonic degrees of freedom dynamically and 
adaptively. The accuracy of the method is first demon¬ 
strated by comparison with reliable results for a nonequi¬ 
librium polaron problem [20]. Then its performance is 
illustrated with a study of wave packet scattering by a 
small EP-coupled structure. 

In quantum lattice models phonon degrees of freedom 
are represented by bosonic sites. As the Hilbert space of a 
single bosonic site is already infinite, it must be truncated 
to a subspace of dimension M from the start in wave- 
function-based numerical approaches [9] . The most com¬ 


mon choice is to use the lowest M eigenstates of a (well 
chosen) boson number operator b^b defining a bare boson 
basis. Then exact diagonalization or exact time propa¬ 
gation can be easily performed but the computational 
cost increases very rapidly with M and exponentially 
with the number of lattice sites. Matrix-product-state 
(MPS) algorithms, such as the density-matrix renormal¬ 
ization group (DMRG) [21-25] and TEBD, allow us to 
treat ID systems at a lower computational cost and thus 
to investigate much larger systems. However, the com¬ 
putational effort still increases as M 3 . Therefore, most 
applications have been restricted to problems with small 
phonon fluctuations (M < 10), in particular for nonequi¬ 
librium problems [26]. 

Instead of a bare boson basis of dimension M. one can 
describe a quantum state \ip) using an optimal local basis 
of dimension Mo < M, which is defined as the eigenba- 
sis of the reduced density matrix of |t/>) for the bosonic 
site [19]. This approach is very efficient for ground-state 
calculations because a sufficient accuracy can be reached 
with a small optimal basis even when a very large bare 
basis would be required. As the optimal basis must be 
calculated self-consistently, the total computational cost 
still rises with M but only linearly. Thus ground-state 
calculations can be carried out with exact diagonaliza¬ 
tion or MPS algorithms for systems with M > 10 3 using 
only moderate computer resources [9, 19, 27, 28]. 

The LBO has never been combined successfully with 
MPS methods to study EP systems out of equilibrium 
but for a very recent study of the spin-boson model [29] . 
The key problem is that the local optimal basis depends 
on the represented quantum state | ip(t)) and thus evolves 
with time. This is clearly seen in our recent study of the 
optimal boson basis for a nonequilibrium polaron prob¬ 
lem (see Figs. 18-20 in Ref. [20]). Therefore, we have de¬ 
veloped an algorithm which allows us to optimize the lo¬ 
cal basis dynamically for the evolving target state \ip{t)). 

We have implemented this approach within the TEBD 
algorithm [17, 18] which is one of the simplest time- 
dependent MPS methods [30-32], For a chain with L 
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FIG. 1. (Color online) Graphical representation of the MPS 
(a) in the original TEBD algorithm and (b) in the TEBD- 
LBO algorithm. 


sites the MPS representation of a quantum state |t/>) in 
an occupation number basis is 

V>(fci,...,fci) = r 1 ,fel A 1 r 2 - fc 2 A 2 (i) 

where the indices kj label the basis states of the dj 
dimensional Hilbert space representing the degrees of 
freedom on the lattice site j £ {1,...,L}. (For a 
bosonic site, dj = M.) The entanglement between two 
parts of the lattice (e.g., the sites {1,..., j} and the sites 
{j + 1 is encoded in the Dj -dimensional posi¬ 

tive definite diagonal matrices A 7 . Hence the matrices 
have dimensions Dj -1 x Dj (with D$ = Dl = 1). 
This MPS is represented graphically in Fig. 1(a). We 
call D = max {D i,..., D^_i} the bond dimension of the 
MPS and d = max {d \,..., dz,} its local dimension. 

Using orthogonality relations for the matrices T and A, 
the matrix elements of the reduced density matrix p 7 for 
the site j are given by 


modifies the two matrices T and the one matrix A as¬ 
sociated with a bond. In practice, we use a second-order 
TSD resulting in an error 0{t 3 ) per time step or 0(t 2 ) 
for a finite period of time. 

Here we only explain how we perform the local update 
in our TEBD-LBO algorithm as our method is otherwise 
identical to the original TEBD [17, 18]. We assume that 
we know the MPS (1) at a given time t in its truncated 
optimal local basis, i.e., we know the To, A and R ma¬ 
trices with dj < do < d. A local update for a single 
bond consists of four steps. First, we build the rank-four 
tensor 


= [\ j - 1 r^\irj+ 1 ’ k *+'v+ 1 ] ai!a2 , (4) 

where the T matrices are given by (3) and the indices 
a 1,2 = 1,..., D number the matrix rows and columns. 
Then, in a second step we carry out the time evolution 
as done in the original TEBD algorithm, 
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( 5 ) 


where U denotes the djdj+i x djdj+± matrix represen¬ 
tation of the local time-evolution operator e~ lTHj ’ j+1 in 
the bare basis. Generally, the computational cost for this 
step is 0(d' l D 2 ) but it can be reduced to 0(d 2 D 2 ) using 
the sparseness of the matrix representation of (7) in a 
bare boson basis. In the third step, we compute the local 
reduced density matrix (2) from the tensor cp using the 
relation 



(A - 7-1 ) 2 (A 7 ) 2 


( 2 ) 


The eigenbasis of this dj x dj matrix is called the optimal 
local basis. The unitary transformation from the optimal 
to the bare basis representation is denoted by i? 7 


do 

= E R LsA Si - 

So —1 


( 3 ) 


This transformation is exact if do = dj. The matrices 
A- 7 are not affected by the basis change. The new MPS 
structure is illustrated in Fig. 1(b). This transformation 
of the matrices T in the TEBD algorithm is similar to 
the approach proposed for a variational MPS [28, 33]. 
Each optimal basis state has a weight (eigenvalue) in 
the interval [0,1]. We can thus approximate the origi¬ 
nal state (1) using only the do(< d) eigenstates with the 
highest weights. 

For a Hamiltonian which includes only on-site and 
nearest-neighbor interactions H = y~h JJ,.,_n, the time- 
evolution of the MPS (1) can be decomposed into suc¬ 
cessive local updates with a time step r using a Trotter- 
Suzuki decomposition (TSD) of the time evolution oper¬ 
ator: ip = e~ lTHj - j+1 ip, where the local operator H ;l j +1 
acts only on a single bond (i.e., sites j and j +1) [17, 18] . 
Each local update is a unitary transformation which 
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and then diagonalize it. This yields the new optimal 
bases for the sites j and j + 1, i.e., new transformations 
i?- 7 and R 1+l . We discard the optimal eigenstates with a 
negligible weight (e.g., lower than 10 -13 ) and thus obtain 
a new truncated optimal basis of dimension do < d. The 
tensor cp is then projected into the new optimal basis. 
Finally, in a fourth step the new matrices f^,, Pq” 1 , and 
A- 7 are calculated from the projected tensor (po exactly 
as in the original TEBD algorithm. Hence we obtain the 
MPS representation of the state ip at time t + r in its 
optimal local basis. 

In summary, we repeatedly propagate the wave func¬ 
tion (1) in a bare local basis to enlarge the effective 
Hilbert space and then project it onto a new effective 
Hilbert space using the LBO to control the dimension of 
the MPS. If do < d, the total computational effort scales 
as d 3 D 3 exactly as with a bare basis. If a small optimal 
basis is sufficient (do ^ d ), however, our algorithm scales 
as the largest of d 3 Q D 3 (the computational cost of TEBD 
in the optimal basis) and d 3 D 2 (the computational cost 
for adapting the optimal basis) and thus it is significantly 
faster than a bare basis simulation. Contrary to the lin¬ 
ear scaling of ground-state methods with d [9, 19, 28], 
however, the computational cost still increases as d 3 but 
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the prefactor is reduced significantly, in particular by a 
factor oc 1/D. Therefore, the advantage of the TEBD- 
LBO will become more pronounced for problems with a 
large block entanglement, i.e., D 1. 

Next we turn to two applications of our TEBD-LBO 
method to EP systems out of equilibrium. We consider 
an Lh ~site Holstein chain [ 8 ] connected at each end to 
tight-binding leads with Ltb fermion sites. The Hamil¬ 
tonian of the full system (with L = Lh + 2Ltb sites) 
is 


L — l 

H = ~ t oJ2 ( c k+i+ c 3+i c i) ( ? ) 

3 =1 

Ltb+Lh 

+ £ 

3 = I j tb 1 

where bj and Cj annihilate a phonon (boson) and a (spin- 
less) fermion on site j, respectively, and rij = cjc^. Thus 
d = 2 M in this model. The model parameters are the 
phonon frequency wo > 0, the EP coupling 7 and the 
hopping integral to- We work with h = 1 and set the 
energy scale to = 1 , thus the time unit is h/to = 1 . 

Here we restrict ourselves to the nonequilibrium dy¬ 
namics of an electron coupled to phonons (i.e., the po- 
laron dynamics), which has recently become a widely 
studied topic [20, 34-45]. One-electron problems have 
MPS with low bond dimensions D , which can easily be 
simulated on a workstation when the effective local di¬ 
mension is small (d or do 10). Thus they provide us 
with a practical test field for our TEBD-LBO method. 
Typically, we use D < 30 or kept all block eigenstates 
with weight > 10 -15 in combination with a time step r 
as small as 10 -3 to keep TEBD errors (induced by the 
TSD and the truncation of the bond dimensions) under 
control. The conservation of the electron number is used 
to decompose the matrices P, A, p, and <j> into block sub¬ 
matrices and thus speed up the calculations. Therefore, 
we also obtain different optimal boson states as a func¬ 
tion of the electronic occupation of a site. 

The initial wave function contains no phonon 

L 

l^ = O))=^/(j)ct|0), (8) 

j =1 

with the vacuum state |0). Thus it is only slightly en¬ 
tangled, as D = d = 2. Naturally, these dimensions 
can increase significantly when the wave function evolves 
with time [25]. Consequently, we always start our simu¬ 
lations with a small bare basis dimension d. After every 
time step r, d is increased if the occupation of the highest 
phonon state exceeds some threshold (e.g., 10 -7 ). 

First, we test our algorithm on the dynamics of a 
highly excited electron coupled to phonons [20]. In 
that case, no tight-binding lead is attached to the Hol¬ 
stein chain (L = Lh) and the electron is initially in 
a state /(j) = yj2/(L + l)sin(Aj) with, e.g., K = 


bjb j -'y(b] + bAn j 



t 


t 


FIG. 2. (Color online) Comparison of LFS and TEBD-LBO 
for various do- (a) Time evolution of the phonon number 
calculated for L = Lh = 6 , 7 = 2 and ljq = 1. (b) Relative 
deviations of the TEBD-LBO data from the LFS results. De¬ 
viations for TEBD with a bare basis of dimension d = 86 are 
also shown. 


7 rL/(L + 1). We have recently investigated this prob¬ 
lem [ 20 ] using very accurate simulations within a limited 
functional space (LFS) [38, 46] and we use these results 
to check the method presented here. As an illustration, 
Fig. 2 compares the time evolution of the phonon num¬ 
ber Aph = (J2j bjbj) calculated using TEBD-LBO and 
LFS. In Fig. 2(a) we see that the agreement is very good 
while in Fig. 2(b) we observe that relative deviations be¬ 
come overall smaller when the dimension do is increased 
and approach the values obtained using TEBD with a 
bare basis. Therefore, the global exponential increase of 
deviations with time is not due to the LBO but to in¬ 
trinsic numerical errors of the TEBD and LFS methods. 
However, our tests also confirm that for this problem do 
must be a substantial fraction of the bare basis dimension 
d (e.g., do ~ d/4) to achieve a similar accuracy. Conse¬ 
quently, the dynamical LBO does not reduce the compu¬ 
tational cost significantly in comparison to the bare basis 
approach for this type of problem. This is due to the 
relatively broad distribution of the local density-matrix 
eigenvalues that we found in our previous work [ 20 ]. 

Second, we apply our method to the scattering of an 
electronic wave packet by a small EP-coupled structure. 
In that case the tight-binding leads are much longer than 
the Holstein chain (we use Ltb up to 280 and Lh < 6 
sites). The initial state is a Gaussian wave packet cen¬ 
tered around a site jo in the left lead and with a positive 
velocity vq ~ 2 f 0 sin(AT) 


f(j) = Cex p 


(j -jo) 2 ' 

4 a 2 


exp[iKj], 


(9) 


where C is a normalization constant, Ltb — jo 3> a ^ 1, 
and 7 r > K > 0. For the calculations presented here, we 
use cr = 5 and K = 7 r/ 2 . After a time t « (Ltb — jo)/v 
the wave packet reaches the Holstein chain where it be¬ 
comes temporarily self-trapped, and finally it is partially 
transmitted and reflected [47]. 

For this problem we find that the dynamical LBO 
reduces the computational effort substantially when 
bosonic fluctuations are large. For cases which can be 
simulated with both a bare basis and an optimized one, 
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FIG. 3. (Color online) (a) Electronic density ( rij) as a func¬ 
tion of site j and time t for Lh = 6, 7 = 1, and wo = 2.25. The 
Holstein chain corresponds to sites j = 280,..., 285. (b) To¬ 
tal electronic density tlh in a Holstein chain of length Lh = 6 
(red solid line), 3 (blue dotted line) and 1 (black dashed line) 
as a function of time t. The inset shows the same data on a 
logarithmic scale. 


we already observe speed-up factors larger than 10. For 
instance, for wo = 0.6 and 7 = 2 calculations with a 
bare basis of dimension d = 124 (M = 62) take 14 
times longer than with an optimized basis of dimen¬ 
sion do = 9, but both approaches yield similar results 
with relative deviations smaller than 10~ 3 . For larger 
phonon fluctuations, we can only complete simulations 
using TEBD-LBO. For instance, in the strong-coupling 
adiabatic regime (uq = 0.2 and 7 = 2), the required bare 
basis dimension is of the order of 10 3 but we can perform 
the TEBD-LBO simulation using only up to do = 23 op¬ 
timal states. Therefore, the dynamical LBO allows us to 
study regimes that we could not treat with the standard 
TEBD algorithm on our workstation. (For comparison, 
dimensions d = 30 and D = 5 were reported in Ref. [29].) 

The direct injection of an electronic wave packet into 
an EP-coupled chain was studied previously [34, 37] but 
the scenario considered here has not been studied yet. 
Thus we briefly discuss two interesting phenomena that 
we have observed but postpone a more thorough discus¬ 
sion to a future work. The first phenomenon is the tem¬ 
porary self-trapping of the electron in the EP-coupled 
structure. In Fig. 3(a) we see that the electron reaches 
the EP-coupled structure at t « 30 but that a finite den¬ 
sity remains in that region even for t = 150 [47]. At 
several times, fractions of the wave packet leave the Hol¬ 
stein chain and start to propagate in the leads. The 
probability of finding the electron in the Holstein chain 

uh = X[^=LTjf+i( n l) ' s shown in Fig. 3(b). It increases 
rapidly when the wave packet reaches the left edge site 
of the Holstein chain at t s=s 30 and then decays expo¬ 
nentially fast for longer times. The decay rate is longer 
for longer chains. Therefore, (a fraction of) the elec¬ 
tronic wave packet becomes temporarily self-trapped in 
the EP-coupled structure and is belatedly transmitted or 
reflected. 

The second phenomenon is the dissipation of the elec¬ 
tron energy due to inelastic scattering processes. Fig¬ 
ure 4(a) shows that one pair of transmitted and re¬ 
flected wave packets moves with the same absolute ve¬ 
locity vo = 2 as the incident wave packet while a sec¬ 




FIG. 4. (Color online) (a) Electronic density ( rij } as a func¬ 
tion of site j and time t for Lh = 1,7 = 2.5, and u>o = 1.65. 
The position of the EP site is j = 90. (b) Phonon energy 
E p h = woNph as a function of time for Lh = 1, wo = 1.5, and 
several values of 7. 

ond pair moves with a lower velocity v\ m 1.1 [47]. 
This corresponds to an inelastic process where a phonon 
is excited by the presence of the electron in the EP- 
coupled structure and then left behind when the electron 
propagates away from this structure. The final veloc¬ 
ity v n = 2fosin(fc rl ) can easily be determined from the 
equality of the asymptotic total energy for t —> ±00, 
— 2to cos (K) = ncoo — 2to cos (k n ), where n is the number 
of excited phonons left behind. Similar patterns have 
been observed recently in a ID photonic wave guide cou¬ 
pled to a two-level scatterer, see Fig. 3 in Ref. [48]. In 
Fig. 4(b) we see that the phonon energy E p h = wofVph, 
which is zero initially, remains finite after the electron 
has left the EP-coupled structure (i.e., for t —>• 00). This 
confirms that an irreversible energy transfer occurs from 
the electron to the phonon degrees of freedom (dissipa¬ 
tion) . 

In summary, we have developed a TEBD algorithm 
with a dynamical optimization of local boson bases that 
allows us to simulate the nonequilibrium dynamics of 
electron-phonon systems more efficiently. This opens 
the way for numerous theoretical investigations of time- 
resolved spectral properties [1, 2], photoinduced phase 
transitions [3, 4], and transport in low dimensions [5-7]. 
The overall performance depends on the properties of the 
local density matrix out of equilibrium and thus on the 
specific problem investigated. The basic idea can eas¬ 
ily be combined with other time-dependent MPS meth¬ 
ods [30-32] or applied to other bosonic systems such as 
correlated photons [12, 13], quantum baths [14, 15], and 
scalar fields [16] as well as to other systems with large 
local Hilbert spaces such as high spin models [49, 50]. 
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